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ABSTRACT 

We present a data simulation package designed to create a series of simulated 
data samples for a detector with non-destructive sampling capability. The original 
(yr) ' intent of this software was to provide a method for generating simulated images for 

the Next Generation Space Telescope, but it is general enough to apply to almost 
I/-) ■ any non-destructive detector or instrument. MultiDataSim can be used to generate 

"practice" sampling strategies for an instrument or field-of-view, and thus can be used 
. to identify optimal observation strategies. 

O 1. Introduction 

MultiDataSim was written to test a cosmic-ray identification-and-rejection program 



(Offenberg, et. al., 1999), written for the Next Generation Space Telescope (NGST). However, it 
is useful for many purposes and for other observatories. One use of MultiDataSim is to compare 



on the merits of various data sampling techniques and timing approaches when presented with 
a detector capable of non-destructive reads. By reading the detector many times, it is possible 
to manipulate the data in ways that would otherwise be unavailable, but defining an optimum 
strategy can become complicated. A study of this trade-space (Fixsen, et. al., in press) suggests 
that there is no general solution to this question: it depends on the observing program and the 
specific parameters of the observation. 

This paper is divided into three parts. The first is a narrative background description of 
MultiDataSim and its potential uses. The second section provides the documentation for the 
software. The final section, the appendix, contains the C++ source code for the software itself. 
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These studies are supported by the NASA Remote Exploration and Experimentation Project 
(REE), which is administered at the Jet Propulsion Laboratory under Dr. Robert Ferraro, Project 
Manager. 

2. Background 

There are several methods for observing a single field-of-view when a non-destructive 
capable-detector is used. We discuss two of them here (more information can be found in Fixsen, 
et. al., in press). 

The first is Fowler Sampling (Fowler & Gatley, 1990), in which the detector is sampled N/2 
times at the start of the observation, followed by a long integration and is then sampled N/2 times 
at the end of the observation. Essentially, we are measuring the accumulating charge N/2 times, 
so we can reduce the random read noise introduced by the detector by a factor of y/N/A. 

A second method is to sample the data Up-the-Ramp, at set intervals during the observation — 
sampling at uniform intervals being the simplest and most straightforward case. The signal is 
then computed as the slope of the ramp via a least-squares fit. The read-noise reduced by a factor 
of y/N/\2 (in the uniform- weighted case). Up-the-Ramp sampling also allows observers to search 
the ramp for glitches (e.g. those caused by cosmic rays), saturated pixels, drift in the tracking, 
transient jitter and so on. Such suspect data could then be excluded from the processing while 
preserving the rest of the data. As this is just one possible trade, it is easy to see that Fowler vs 
Up-the-Ramp vs some other sampling method can be a complicated trade space. 

The MultiDataSim tool allows the user to select an arbitrary set of sampling times; Uniform 
and Fowler sampling algorithms are already built into the code. The software starts with a 
noiseless "Real Sky" image and generates simulated observations for the specific intervals with 
appropriate noise contamination. Each sample is affected by Poisson noise (i.e. photon-counting 
statistics), random read noise and cosmic rays. Most of the default values are chosen for the 
NGST. Thus, the cosmic ray rate (5 particles cm~ 2 s _1 ) is that expected for a deep-space (L2) 
orbit, the physical pixel size is the baseline for NGST (27^ x 27 '\x x 10/z), and so on. It is possible 
to turn off these features, and override the default values on the command-line. The output images 
are generated as a sequence of FITS (Flexible Image Transport System) files, one sample of the 
detector array per file. 

There are two functions which the authors have not yet had the opportunity to implement in 
MultiDataSim, although work-arounds exist for each. The first is to apply a non-linearity curve 
to the data, although this can be applied after the fact. MultiDataSim also does not apply a 
point-spread function (PSF) to the data. This can easily be done by applying the PSF to the 
input image. 
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3. Documentation 

Calling sequence: 

MultiDataSim [options] inputfile [options] 

where inputfile is a FITS file containing the "ideal" input image. Most of the options are 
order-independent and can appear before or after the inputfile argument. 

Options: 

-co x Specifies a constant DC offset (versus dark current) to be applied to the 

image (default: 0.0). 

-craa Specify to use cosmic-ray area-of-effect. This turns on a small mask so 1% 

of the cosmic ray flux ends up in each of the 4 adjacent -neighbor pixels and 
.1% into the 4 diagonal-neighbor pixels. 

-err x Specifies cosmic ray rate as x events cm -2 sec -1 (default: 4). 

-de x Specifies the DC dark current is x e~ sec -1 (default: 0.02). 

-df file Specifies the path to the image containing the dark frame (default: dark 

frame is blank). 

-expf file Specifies that the software should use a customized sample times instead of 

either Fowler or Uniform sampling, file is the name of an ascii file containing 
a list of sample times (one per line). The software will record a t=0 image 
as well as the ones specified. Using the -expf keyword overrides the "-fs", 
"-n", and "-t" keywords. 

-ff file Specify Flat field image (default: none.) 

-g x Specify the gain in e~ data-unit -1 (default: 4). 

-fs Specifies that the software should use Fowler Sampling instead of the default 

Uniform Sampling method. The interval between Fowler samples is 5 
seconds. 

-ie x Specifies the exposure time of the "ideal" input image in seconds (default: 

1000). 

-n x Specifies that the observing sequence contains x observations including an 

observation at t=0 seconds (default: 64). 

-nc Specifies that the software should not inject cosmic rays into the detector. 

-ncf Specifies that the software should not record the injected cosmic ray events 

in separate files (see "-o"). 

-np Specifies that the software should not apply Poisson noise to the data. 
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-nr Specifies that the software should not apply readnoise to the data. 

-o file Specifies the base file name for the output sequence (default: "datasim_"). 

Each simulated data sample is stored in its own FITS file named 
fileXXXXXX.fit, where XXXXXX is the time of the sample during the 
integration, in whole seconds. Note that this notation imposes 2 limits: 
elapsed time between samples must be a minimum of 1 second, and the 
observing sequence as a whole is limited to 10 6 seconds. In addition to 
creating fileXXXXXX.fit, datasim also records cosmic ray events in a series 
of files fileXXXXXX.fit_cosmic (unless the -ncf argument is specified). 



-pix l,w,d Specify length, width, depth of pixel in fi (default: 27,27,10). 

-r x Specify read noise in e~ (default: 15). 

-sc x Specifies angular shield coverage is x degrees (default: 1.0). 

-t x Specifies that the observing sequence lasts x seconds (default: 1000). 

-v Specifies that the software should run in "Verbose" mode. 
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APPENDIX 

The code here is long enough as it is, so we have omitted a few functions and files from this 
appendix. We include the header files which define C++ classes, but omit all other header files. 
The other header files contain just procedure prototype lines, so recreating them is straightforward. 

The rani function (called in cosmic. cc) and the poidev and gasdev functions (called in 
detectorimage.ee) are all from the Numerical Recipes library, and can be found in the text 
Numerical Recipes (Press, Flannery, Teukolsky & Vetterling). We are using a version of the 
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on-line library based on the 1986 edition of Numerical Recipes; unfortunately, that version 
of the library is Fortran-based, so the translation to C was awkward in places. We created a C 
library, librecipes_g+- 1-, containing rani, poidev and gasdev, plus the Numerical Recipes functions 
upon which these 3 depend. In addition to the Numerical Recipes library, we are also using the 
standard C math library. 

Implementing FITS I/O is left as an exercise for the reader. The -lefitsio in the Makefile 
refers to the "CFITSIO" library. The CFITSIO library is available from other sources on the web 
and is not included in this paper. We are using version 1.4, by William Pence (1998). 

A. Makefile 

# MAKEFILE for MultiDataSim 
# 

# Makefile to generate the MultiDataSim software. This Makefile was 

# written for the standard Sun Solaris (tm) Make utility. 
# 

# Written by Joel D. Offenberg, Raytheon ITSS 
# 

OBJ = datasim.o det ect or image . o fits.o string2num.o cosmic. o 

INC_DIR = -I. ./include 

LIB_DIR = -L../lib -L/usr/local/lib 

LIBS = -lefitsio -lm -lrecipes_g++ 

FLAGS = -0 

CCC = g++ $ (FLAGS) 

MultiDataSim : $(0BJ) 

$(CCC) -o MultiDataSim $(0BJ) $(INC_DIR) $(LIB_DIR) $(LIBS) 

install : MultiDataSim 
mv MultiDataSim ../bin 
chmod 755 . . /bin/MultiDataSim 

. cc . o : 

$(CCC) -c -o $*.o $*.cc $(INC_DIR) 
# 
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B. pixel.h 

/ ***** 

PIXEL.H 

This is a simple class to encapsulate several properties of a pixel in 
the detector. At this writing, it is just the size, although we might 
end up encapsulating other properties over time for neatness' sake. 

Written by Joel Dffenberg, Raytheon ITSS 

******/ 

#ifndef PIXEL_H 
#define PIXEL_H 

class PIXEL 
{ 

public : 
float x, y, z; 

PIXEL () { x = 27.; y=27.; z=10.; } 

PIXEL(float a, float b, float c) { x = a; y = b; c = z; } 

}; 

#endif 



C. detectorimage.h 

/ ****** 

detectorimage . h 

File containing the class definitions for the Detectorimage class. 

Written by Joel D. Dffenberg, Raytheon ITSS 

*****/ 

#ifndef DETECTOR_IMAGE_H 



#define DETECTOR_IMAGE_H 



#include "pixel. h" 

class Detector Image 
{ 

private : 

float ** image; 

long axisl; 

long axis2; 

float exptime; 

float deltatime; 
public : 

DetectorImage(int , int) ; 

Detectorlmage (float *, int, int, float); 

DetectorImage(DetectorImage* ); 

Detectorlmage (char *f ilename) ; 

Detectorlmage (char *filename, long startrow, long started, 
long numrow, long numcol) ; 
"Detectorlmage () ; 

void zero_out(); 

void setval(int i, int j, float value) { image[i][j] = value 
float getval(int i, int j) { return image[i][j]; } 

float *Buffer() ; 

long Axisl () { return axisl; } 

long Axis2() { return axis2; } 

float exposure () { return exptime; } 

void set_exposure (float exp) { exptime = exp; } 

float get_deltatime () {return deltatime;} 

void scale_to_exptime (float t) ; 

void count ing_st at istics () ; 

void readnoise (float) ; 

void roundoff (float) ; 

void roundoff () { this->roundof f (1 . ) ; } 

void write(char *f ilename) ; 
void write_int (char *f ilename) ; 
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void plus_equals (Detectorlmage) ; 

void plus_equals (Detectorlmage *) ; 

void plus_equals (float) ; 

void divide_equals (float c) ; 

void times_equals (float c) ; 

void times_equals (Detectorlmage *A) ; 

void times_equals (Detectorlmage A) {times_equals (&A) ; } ; 
void ceiling(f loat c) ; 

void add_cosmicrays () 
{ add_cosmicrays (exposure () ) ; } 

void add_cosmicrays (float time) 
{add_cosmicrays (time , 1.0);} 

void add_cosmicrays (float time, float shieldcover) 
{add_cosmicrays (time , shieldcover ,4. ) ; } 

void add_cosmicrays (float time, float shieldcover, float crrate) 
{add_cosmicrays (time , shieldcover , crrate , 0) ; } 

void add_cosmicrays (float time, float shieldcover, float crrate, char mask) 
{add_cosmicrays (time , shieldcover, crrate, mask, *(new PIXEL()));} 

void add_cosmicrays (float time, float shieldcover, float crrate, 
char mask, PIXEL P) ; 

void test_f or_inf inity (char *message) ; 

}; 

#endif 



D. detectorimage.ee 

/ ***** 

detectorlmage . c 

C code files for the Detectorlmage class. See detectorimage .h. 

Written by Joel D. Dffenberg, Raytheon ITSS 
*****/ 

#include <stdio.h> 



-9- 



#include <string.h> 
#include <math.h> 
#include <nr.h> 
#include "fitsio.h" 
#include "pixel. h" 
#include "string2num.h" 
#include "fits.h" 
#include "detectorimage .h" 
#include "cosmic. h" 

extern long idnum; 

/** Constructors **/ 

/* Create an empty image naxisl x naxis2 */ 
Detectorimage :: Detectorimage (int naxisl, int naxis2) { 
int j,i; 

image = new (float*) [naxisl] ; 
for (j=0; j<naxisl; j++) { 

image [j] = new float [naxis2] ; 

for (i=0; i<naxis2; image [j] [i++] = 0.0); 

> 

axisl = naxisl; 
axis2 = naxis2; 
exptime = 1000. ; 

} 

/* Create an image and fill with naxisl x naxis2 elements from ' 'buffer.' J */ 
Detectorimage :: Detectorimage (float *buffer, int naxisl, int naxis2, float exp) { 

int i , j ; 

long count=0; 

image = new (float*) [naxisl] ; 
for (j=0; j<naxisl; j++) { 

image [j] = new float [naxis2] ; 

for (i=0; i<naxis2; image [j] [i++] = buf f er [count++] ) ; 

> 
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axisl = 
axis2 = 
exptime 



naxisl ; 
naxis2 ; 
= exp; 



/* Make a duplicate image */ 

Detector Image :: Detector Image (Detector Image *A) { 
float *buffer; 
int i , j ; 
long count=0; 

image = new (float*) [A->Axisl ()] ; 
for (j=0; j<A->Axisl() ; { 

image[j] = new float [A->Axis2Q] ; 

for (i=0; i<A->Axis2 () ; { 
image [j] [i] = A->image[j] [i] ; 

} 

} 



axisl = A->Axisl(); 
axis2 = A->Axis2(); 
exptime = A->exposure () ; 



/* Read image from a FITS file on the disk */ 
Detector Image :: Detector Image (char *filename) { 

float *buffer; 

long naxis[2], count=0; 

int i , j ; 

char *saveexp; 

if ( (saveexp = readheader (filename , "EXPOSURE" ) ) == NULL) { 
if ((saveexp = readheader (filename , "EXPTIME") ) == NULL) { 

exptime = 1.0; 
} else { 

saveexp = strchr (saveexp, '=') ; 

exptime = string2num (saveexp) ; 

} 

} else { 

saveexp = strchr (saveexp ,'=') ; 
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exptime = string2num(saveexp) ; 

> 

if ((saveexp = readheader (filename , "DELTATIM" ) ) != NULL) { 

saveexp = strchr (saveexp, ' = ' ) ; 

deltatime = string2num (saveexp) ; 
} else deltatime = 1.0; 
buffer = readimage (filename, naxis) ; 
image = new (float*) [naxis [0] ] ; 
for (j=0; j<naxis[0]; j++) { 

image [j] = new float [naxis [1] ] ; 

for (i=0; i<naxis[l]; image [j] [i++] = buffer [count ++]) ; 

> 

axisl = naxis [0]; 
axis2 = naxis [1] ; 

delete buffer; 

} 

/* Read a section from a FITS file on the disk. */ 

Detector Image :: Detector Image (char *filename, long startrow, long startcol, 
long numrow, long numcol) { 

float *buffer; 

long naxis [2], count=0; 

int i , j ; 

char *saveexp; 

if ( (saveexp = readheader (filename , "EXPOSURE") ) 
if ((saveexp = readheader (filename , "EXPTIME") ) 

exptime = 1.0; 
} else { 

saveexp = strchr (saveexp, '=') ; 

exptime = string2num (saveexp) ; 

} 

} else { 

saveexp = strchr (saveexp, '=') ; 
exptime = string2num (saveexp) ; 

} 



== NULL) { 
== NULL) { 
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if ((saveexp = readheader (filename , "DELTATIM" ) ) != NULL) { 

saveexp = strchr (saveexp , ' = ' ) ; 

deltatime = string2num (saveexp) ; 
} else deltatime = 1.0; 

buffer = readsubimage (filename, naxis, startrow, startcol, numrow, numcol) ; 
image = new (float*) [naxis [0] ] ; 
for (j=0; j<naxis[0]; j++) { 

image [j] = new float [naxis [1] ] ; 

for (i=0; i<naxis[l]; image [j] [i++] = buffer [count ++]) ; 

> 

axisl = naxis [0] ; 
axis2 = naxis [1] ; 

delete buffer; 



/*** Destructor ****/ 
Detector Image : : ~DetectorImage () { 
int i ; 

for (i=0; i<this->Axisl () ; i++) { 
delete image [i] ; 

> 

delete image; 

}; 



/*** Buffer ***/ 

/* Returns a pointer to a single buffer which contains the entire image in 
one long stretch instead of the standard 2-D array that the image'' 
contains . 

*/ 

float* Detector Image :: Buffer () { 
float *result; 
int j ; 
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result = new float [axisl*axis2] ; 

for (j=0; j<axisl; j++) { 

memcpyC (result+j*axis2) , image[j], axis2*sizeof (float) ) ; 

} 

return result ; 



/*** plus_equals ***/ 
/* 

Performs the += operation on two Detector Images or a detector image 
plus an offset. I did it this way (instead of overloading the += operator) 
because it removed an ambiguity when dealing with a pointer to a 
Detectorlmage . 

*/ 

/* Plus_equals with an image "a". */ 

void Detectorlmage: :plus_equals (Detectorlmage a) { 

int i , j ; 

if ((a.AxislQ != this->Axisl () ) || (a.Axis2() != this->Axis2 () ) ) { 
printf("\n These images are not the same size!"); 
return; 

} 

for (i=0; i<this->Axisl () ; i++) { 
for (j=0; j<this->Axis2() ; j++) { 
this->image [i] [j] += a. image [i] [j] ; 

} 

} 

exptime += a. exposure () ; 



/* Plus_equals with a reference to an image "*a". */ 
void Detectorlmage: :plus_equals (Detectorlmage *a) { 
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int i,j; 

if ((a->Axisl() != this->Axisl () ) || (a->Axis2() != this->Axis2() ) ) { 
printf("\n These images are not the same size!"); 
return; 

> 

for (i=0; i<this->Axisl () ; i++) { 
for (j=0; j<this->Axis2() ; j++) { 

this->image [i] [j] = this->image [i] [j] + a->image [i] [j] ; 

} 

> 

exptime += a->exposure () ; 

} 

/* Plus_equals with a single value. */ 

void Detectorlmage: :plus_equals (float offset) 

{ 

int i,j; 

for (i=0; i<this->Axisl () ; i++) { 
for (j=0; j<this->Axis2() ; j++) { 
this->image [i] [j] += offset; 

} 

> 

} 

/**** Detectorlmage :: roundoff ****/ 
/* 

Rounds off image data to the nearest n-th integer. Also cuts off the maximum 
value at 2" 16, to prevent 16-bit integer overflow. 
*/ 

void Detectorlmage :: roundoff (float n) { 
int i,j; 

for (i=0; i<this->Axisl () ; i++) { 



-15- 



for (j=0; j<this->Axis2() ; j++) { 

image [i][j] = (image [i] [j] >65535 .) 765535 . : (float) ( (int) image [i] [j] /n)*n; 

} 

> 

} 

/**** Detector Image :: write ****/ 

/* Saves the image to disk as a FITS file. */ 

void Detectorlmage :: write (char *filename) { 
long naxis [2] ; 
float *buffer; 

buffer = this->Buff er() ; 
naxis [0] = this->Axisl () ; 
naxis [1] = this->Axis2 () ; 

writeimage(buf f er , this->exposure () , filename, naxis); 
delete buffer; 

return; 

} 

/**** Detectorlmage : :write_int ****/ 

/* Saves the image to disk as an INTEGER FITS file. */ 

void Detectorlmage : :write_int (char *filename) { 
long naxis [2] ; 
float *buffer; 

buffer = this->Buff er() ; 
naxis [0] = this->Axisl () ; 
naxis [1] = this->Axis2 () ; 

writeimage_int (buffer, this->exposure() , filename, naxis); 
delete buffer; 



return; 

} 



- 16 - 



/**** divide_equals ****/ 

/* A routine to divide the image (and exposure time) by a scalar constant. 
It was done this way (instead of overloading /=) to avoid ambiguity 
with pointer arithmetic. 
*/ 

void Detectorlmage :: divide_equals (float c) { 
int i , j ; 

for (i=0; i<Axisl(); i++) { 
for (j=0; j<Axis2(); j++) { 
image [i] [j] /= c; 

} 

} 

exptime /= c; 

} 

/**** times_equals ****/ 

/* A routine to multiply the image (and exposure time) by a scalar constant. 
It was done this way (instead of overloading *=) to avoid ambiguity 
with pointer arithmetic. 
*/ 

void Detectorlmage: :times_equals (float c) { 
int i , j ; 

for (i=0; i<Axisl(); i++) { 
for (j=0; j<Axis2(); j++) { 
image [i] [j] *= c; 

} 

} 

exptime *= c; 

} 

/**** times_equals ****/ 

/* A routine to multiply the image by another image. 

It was done this way (instead of overloading *=) to avoid ambiguity 

with pointer arithmetic. 

*/ 
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void Detectorlmage: :times_equals(DetectorImage *A) { 
int i,j; 

for (i=0; i<Axisl(); i++) { 
for (j=0; j<Axis2(); { 

imaged] [j] *= A->getval(i, j) ; 

} 

} 

> 

/**** scale_to_exptime ****/ 

/* A routine to scale the image and exposure time to the specified number 
of seconds of exposure time. Basically, this uses divide_equals to 
do the division, but it calculates the factor based on the existing 
exposure time rather than making the user do it. 

*/ 

void Detectorlmage :: scale_to_exptime (float t) { 
t = this->exposure () / t; 
this->divide_equals (t) ; 

> 

Detectorlmage *emptyimage; 

/**** add_cosmicrays (float time, float shieldcover, float rate) ****/ 
/* Routine to add cosmic rays to the image. Asks the user to supply the 
number of seconds of integration for the cosmic rays. Ultimately, 
the proper routine will be add_cosmicrays () , which will extract the 
current exposure time and use that to calculate the number of CRs. 
*/ 

void Detectorlmage :: add_cosmicrays (float time, float shieldcover, 
float crrate, char mask, PIXEL P) { 



float **crimage; 
int i,j; 
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float tmp; 

/* Initialize CRImage */ 

crimage = new (float *) [this->Axisl ()] ; 

for (j=0; j<this->Axisl() ; j++) { 

crimage[j] = new (float) [this->Axis2()] ; 

for (i=0; i<this->Axis2 () ; i++) { 
crimage [j] [i] = 0.0; 

} 

} 



/* Compute cosmic rays. We'll need to tweak some of the instrumental 
parameters for cosmic () ; see cosmic. c. 
*/ 

cosmic (crimage ,this->Axisl () ,this->Axis2() , time, shieldcover,crrate,mask, P) ; 

for (i=0; i<this->Axisl () ; i++) { 
for (j=0; j<this->Axis2() ; j++) { 
image [i] [j] += crimage [i] [j] ; 

emptyimage->setval (i , j , empty image->get val (i , j)+crimage[i] [j] ) ; 

} 

} 

for (j=0; j<this->Axisl() ; j++) { 
delete crimage [j]; 

} 

delete crimage; 



/*** counting_statistics ***/ 

/* Takes "ideal" image and generates a version based on Poisson counting 
statistics. */ 

void Detectorlmage : : counting_statistics() { 
int i,j; 



for (i=0; i<Axisl(); i++) { 
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for (j=0; j<Axis2(); { 

image [i][j] = poidev(image [i] [j] , feidnum) ; 

} 

> 

} 

/*** readnoise ***/ 

/* Takes "ideal" image and generates a version based on Gaussian distribution 
statistics. */ 

void Detectorlmage :: readnoise (float variance) { 
int i,j; 
float tmp; 

for (i=0; KAxisK); i++) { 
for (j=0; j<Axis2(); { 

image[i][j] += (float) ( (int)gasdev(&idnum) * variance); 

} 

} 

} 

/*** zero_out ***/ 

/* Simple routine to erase an image. Also sets exptime to zero. */ 

void Detectorlmage :: zero_out () { 
int i ; 
int j ; 

for (i=0; KAxisK); i++) { 
for (j=0; j<Axis2(); j++) { 
image [i] [j] = 0. ; 

} 

} 

exptime = 0.0; 

} 



/* 

Detectorlmage : ceiling 
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Limit value at a top value (i.e. saturation limit). 

*/ 

void Detectorlmage :: ceiling(float c) { 
int i , j ; 

for (i=0; i<Axisl(); i++) { 
for (j=0; j<Axis2(); { 
if (image [i] [j] >= c) { 
image [i] [j] = c; 
} 

} 

} 

} 

/*** test_f or_inf inity ***/ 
/* 

A routine to test if the image contains any ' 'infinite' ' (or 

impossibly large) values. Debugging routine. 

*/ 

void Detectorlmage :: test_for_inf inity (char *message) { 
int i , j ; 

for (i=0; i<Axisl(); i++) { 
for (j=0; j<Axis2(); j++) { 

if (fabs (image [i] [j] > le+20)) { 
printf("°/s \n" , message) ; 

//Short-circuit out of loop if one is found 
i = Axisl () ; 
j = Axisl () ; 
} 

} 

} 

} 
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E. cosmic. cc 

cosmic . cc 

C Code file containing routines to generate simulated cosmic ray hits 
on detector. The entry point for these routines is cosmic: 

cosmic (detector_array , array_xsize, array_ysize, length) 

detector_array [array_ysize] [array_xsize] is a 2-dimensional array to which 

the cosmic rays will be added. 
array_xsize, array_ysize are the dimensions of the detector_array 
length is the number of seconds over which cosmic rays are to be added. 

ASSUMPTIONS : 

1. All CR's liberate 100 (+/- 10) electrons per 0.1 micron travel. 10% 
of the cosmic rays are He nuclei, and thus have 4x the effect. 

2. The plane of the detector (+/- some extra) is blocked by a shield. 

3. The CR sheild stops 100% of CRs, regardless of energy. 

4. A CR event in the detector has no lasting effect. 
Written by: Joel D. Offenberg, Raytheon ITSS 



#include <sys/types .h> 
#include <time.h> 
#include <math.h> 
#include <stdlib.h> 
#include <nr.h> 
#include "pixel. h" 
#include " cosmic. h" 

extern long idnum; 



/*** 
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Global constants 
***/ 

const float deg2rad = 3.14159265/ 180.; /* Number of radians per degree */ 

/ ***** 
cosmic 

Routine to deliver simulated Cosmic ray hits on the detector. 
INPUTS : 

int **image: 2-dimensional array containing image to which cosmic rays 
are to be added. 

int xs, ys: Size of the image, in pixels, to which the cosmic rays are to 
be added. 

float time: Number of seconds of integration for the image, 
float shieldcover Cutoff angle (in degrees) due to shield 

assume shield is perfect 

OUTPUTS : 

The image array is updated. 
Function returns a "1" upon success. 

*****/ 
extern int VERBOSE; 

int cosmic (float **image, int xs, int ys, float time, float shieldcover, 
float rate, char mask, PIXEL P) { 

float edl = 100.; /* # electrons liberated per micron travel of H */ 

float n_rays; /* # cosmic rays on detector */ 
int i,j,k,m,p,q; 

float xpos, ypos, zpos; /* position of cosmic ray hit on detector */ 
float theta, cosphi , sinphi ; /* direction angles (and cosine thereof) of 

cosmic ray velocity */ 

float pathlength; /* Total path length through detector */ 

int npath; /* number of micron steps needed to pass 

through detector */ 

float dpath; /* Number of microns per step. In theory, 
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this is 1 by definition, except for the last 
step, but instead we seem to be assuming 
even not-quite-micron steps. */ 

float depath; /* Number of electrons liberated per step. */ 

float dx, dy, dz; /* Amount of motion per step */ 

float charge; 

int NHe = 0; 

n_rays = xs * P.x * le-4 * ys * P.y * le-4 * rate * time; 

n_rays = poidev(n_rays , feidnum) ; /* Add poisson statistics. */ 
if (n_rays < 0) { n_rays = 0; } /* Correct for negative number */ 

for (i=0; i<n_rays; i++) { 

/* Generate random numbers between 0. . xs-1 and 0..ys-l. */ 
xpos = ranl(feidnum) * (float) xs; 
ypos = rani (feidnum) * (float) ys; 
zpos = 0.0; 

theta = rani (feidnum) *2*3 . 141592653; 

sinphi = sqrt (rani (feidnum) ) ; 
cosphi = sqrt(1.0 - sinphi*sinphi) ; 

if (asin(sinphi)/deg2rad > shieldcover) { 

/* 10% of the time, the CR will be a He nucleus, with 4 times the effect */ 
charge = ( ranl(feidnum) > 0.9) ? 2. : 1.; 
if (charge ==2.) {NHe++;> 

depath = 100. ; 
dpath =0.1; 

dx = cos (theta) * dpath * sinphi; /* x element per step */ 
dy = sin (theta) * dpath * sinphi; /* y element per step */ 
dz = dpath * cosphi; /* z element per step */ 

/* Compute position of cosmic ray at each step and add appropriate 
signal to detector array. Stop when outside the detector. */ 
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for (zpos =0.0; (zpos <= P.z) kk (xpos>=0) kk (ypos>=0) kk 
(xpos<xs) kk (ypos<ys) ; zpos += dz) { 

if (mask == 1) { 
// Add cosmic ray area-of-ef f ect mask. The values are based on 
// Bernie Rauscher's cosmic ray measurements, 
int xxpp,yypp; 

float crmask[3] [3] = {{1 .06e-3, 1 .66e-2, 1 .06e-3>, 

{1.66e-2,1.00,1.66e-2}, 

{1 . 06e-3 , 1 . 66e-2 , 1 . 06e-3}} ; 
float crcont = depath* (charge* charge) ; 

for (xxpp=0; xxpp<=2; xxpp++) { 
for (yypp=0; yypp<=2; yypp++) { 

if ( (xxpp+xpos-1 >= 0) kk (yypp + ypos-1 >= 0) kk 
(xxpp+xpos-1 < xs) kk (yypp + ypos-1 < ys)) { 

image[(int) xpos + xxpp-1] [(int) ypos+yypp-1] += 
poidev(crcont , &idnum) * crmask[xxpp] [yypp] ; 
} > > } 

else 

image[(int) xpos] [(int) ypos] += poidev (depath* (charge* charge) , feidnum) ; 
xpos += dx/P.x; 
ypos += dy/P.y; 

> 

/* If the CR hits the bottom of the detector the odds are it won't hit it 
on an even increment in the steps above. Thus, add an extra term for 
the last little bit. The CR-area PSF is not applied here, since this 
is just a little bit at the end. [It probably should be, but we'll 
live for now.] 

*/ 

if ((xpos>=0) kk (ypos>=0) kk (xpos<xs) kk (ypos<ys)) { 
image[(int) xpos] [(int) ypos] += 
poidev (depath* ( charge* charge ) * (P. z-zpos+dz)/dz, feidnum) ; 
} 
} 

> 
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if (VERBOSE) printf ("\t%d CRs, %d He\n" , (int) n_rays, NHe) ; 

/* DONE */ 
return 1 ; 



F. datasim.cc 



/ ***** 
Main. c 



Main file for data NGST warts-and-all data simulator. 



Written by: Joel D. Offenberg, Raytheon ITSS 



*****/ 



#include <string.h> 
#include <sys/types .h> 
#include <time.h> 
#include <stdio.h> 
#include "pixel. h" 
#include "fits.h" 
#include " det ect or image .h' 
#include "string2num.h" 
#include "datasim.h" 



long idnum; 
int VERB0SE=0; 



extern Detector Image *emptyimage; 

/*** main ***/ 
/* 

the entry point and controlling routine. 
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*/ 



int main(int argc, char **argv) { 

char filename [80] ; //input file name 



// Image objects 

Detectorlmage *ldeallmage , *DeltaImage, * Integrate Image , *IntegratePrime , 
*Flat Image ; 



//Other file names. 

char outfilename[80] , outf ilebase [70] , darkf ile [80] ,f latf ile [80] ; 

//exposure time in seconds. 



float exptime; 
float numimg; 
float duration; 
float *deltatimes; 
float offset =0.0; 
float dccurrent=0.02; 
float initexp = 1000; 
float shieldcover = 1 
float crrate = 4. ; 
float rnvar = 15.; 
float gain = 4 . ; 



0; 



//number of intervals, 
//keep track of exposure time to date. 
//List of exposure-time intervals. 
//DC offset 
//Dark current (e-/s) 
//Exposure-time of original image. 
//Angular coverage of shield around plane . 
//Number of cosmic rays/s/cm/cm. 
/ /Readnoise . 
//gain (e-/data unit) 



char cosmicbool=0, readnoise=0, cosmicout=0, exptimef ile [100] ; 

char tempstring [80] ; 

int count ; 

FILE *inputfile; 

PIXEL *P; 



P = new PIXEL () ; 



if (parse_args (argc , argv, filename, outf ilebase, darkf ile, f latf ile, 
feoffset, feduration, fenumimg, fecosmicbool, fereadnoise, fecosmicout, 
fedccurrent, feinitexp, feshieldcover , exptimefile, fecrrate, &rnvar, 
&gain, P)) { 
return 0; 

> 



/* Read in source image. Assume that exposure time is 1000 seconds. */ 
if (strcmp(f ilename , " " ) == 0) { 
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Ideallmage = new Detectorlmage(1024, 1024) ; 
} else { 

Ideallmage = new Detector Image (filename) ; 
IdealImage->set_exposure(initexp) ; 

> 

/* Read in flatfield image. If not specified, make it a constant image. */ 
if (strcmp(flatfile,"") == 0) { 

Flatlmage = new Detectorlmage (IdealImage->Axisl () , IdealImage->Axis2()) ; 

Flatlmage->plus_equals (1 . ) ; 
} else { 

Flatlmage = new Detectorlmage (flatfile) ; 

> 

emptyimage = new Detectorlmage (IdealImage->Axisl () , IdealImage->Axis2() ) ; 

//Seed the random number generator: use the system clock, 
idnum = -time(0) ; 

/* Initialize and zero-out Int egr ate Image . This is a little wasteful 
from a computation standpoint, but it guarantees that Integratelmage 
will be the right size. 

Modified to include possibility that a non-zero dark and/or a constant 
offset should be applied to the image. */ 

if (strcmp (darkf ile, "") == 0) { 

Integratelmage = new Detectorlmage (Ideallmage) ; 

IntegrateImage->zero_out () ; 
} else { 

Integratelmage = new Detectorlmage (darkf ile) ; 

} 

Integratelmage->plus_equals (offset) ; 
//Write out darkf ield image. 

IntegratePrime = new Detectorlmage (Integratelmage) ; 
IntegratePrime->times_equals (Flatlmage) ; 
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if (readnoise == 0) { IntegratePrime->readnoise(rnvar) ; } 
//Modify for gain 

IntegratePrime->divide_equals (gain) ; 
IntegratePrime->set_exposure (0.0); 
IntegratePrime->roundof f () ; 

sprintf (outf ilename , "7 s7 06d.f it" , outf ilebase , (int) 0) ; 
IntegratePrime->write_int (outf ilename) ; 

//Subtract this image from the number of images, "numimg." 
numimg -= 1 ; 

if (cosmicout == 0) { 

sprintf (outf ilename , "7 s7 06d.f it_cosmic" , outf ilebase, (int) 0) ; 
emptyimage->write_int (outf ilename) ; 

} 

/* Figure out deltatimes. If "-fs" (Fowler Sample) option was selected, 
first numimg/2 exposures are 5 second each, followed by 1 big exposure, 
then numimg/2 5 second exposures. Otherwise, the delta times are 
all evenly spaced. */ 

deltatimes = new (float) [numimg] ; 
if ((cosmicbool & 8) != 0) { 

delete deltatimes; 

deltatimes = new (float) [1000] ; 

inputfile = fopen(exptimef ile, "r") ; 

for (count=0; ! feof (inputfile) ; count++) { 

f scanf (inputfile, "7os\n" ,tempstring) ; 

deltatimes [count] = string2num(tempstring) ; 

} 

numimg = count ; 
duration = deltatimes [count] ; 
for (count=( int) numimg; count > 0; 
deltatimes [count — ] -= deltatimes [count-1] ) ; 
f close (inputfile) ; 
} else 

{ if ((cosmicbool & 2) != 0) 
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{ 

for (count = 0; count < (numimg)/2; deltatimes [count ++] = 5.0); 
deltatimes [count++] = duration - 5 . 0* (numimg-1) ; 
for ( ; count < numimg; deltatimes [count++] = 5.0); 

} 

else 
{ 

for (count=0; count<numimg; deltatimes [count ++] = duration/numimg) ; 

> 



/*** Repeat loop for equally-spaced observations through duration. ***/ 

for (count=0,exptime=deltatimes [count] ; count<numimg; 
exptime+= deltatimes [++count] ) { 

if (VERBOSE) printf ("%d seconds\n" , (int) exptime) ; 

/* Calculate ideal image scaled to exposure exptime. */ 
Deltalmage = new Detectorlmage (Ideallmage) ; 
DeltaImage->scale_to_exptime (deltatimes [count] ) ; 
Deltalmage->plus_equals (dccurrent*deltatimes [count] ) ; 

/* Modify image by Poisson counting statistics as detected */ 
if ((cosmicbool & 4) != 4) { 

DeltaImage->counting_statistics () ; 

> 

/* Add cosmic rays */ 
if ((cosmicbool & 1) != 1) { 
emptyimage->zero_out () ; 

DeltaImage->add_cosmicrays (deltatimes [count] , shieldcover , crrate , 
((cosmicbool & 16) != 0)?1:0,*P); 

sprintf (tempstring, "°/ 06d" , (int) exptime) ; 
Deltalmage->test_f or_inf inity (tempstring) ; 

} 

Deltalmage->times_equals (Flatlmage) ; 
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/* Remember the history of the observation until now */ 
Integratelmage->plus_equals (Deltalmage) ; 

/* Copy full observation, add read noise (gaussian distribution, 
variance of 15 electron-units) . 

*/ 

IntegratePrime = new Detectorlmage (Integratelmage) ; 

if (readnoise == 0) { 

IntegratePrime->readnoise (rnvar) ; 

} 

/* Generate filename and save the the observation as it has been 
simulated until now. */ 

IntegratePrime->divide_equals (gain) ; 

IntegratePrime->set_exposure( IntegratePrime->exposure () *gain) ; 
//This is needed because the previous step messes up the exptime. 

IntegratePrime->roundof f () ; 

sprintf (outf ilename , "°/ s°/o06d. f it" , outf ilebase , (int) exptime) ; 
IntegratePrime->write_int (outf ilename) ; 

if (cosmicout == 0) { 

sprintf (outf ilename , "%s / 06d.f it_cosmic" , outf ilebase , (int) exptime) ; 
emptyimage->write_int (outf ilename) ; 

} 

/* Free up allocated space before next loop. */ 
delete Deltalmage; 
delete IntegratePrime; 

} 

return 1; //Done. 

} 

/*** parse_args ***/ 
/* 
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Collect all of the argument-parsing routines into one place. 
There are C-library routines that will do most of this, but, 
unfortunately, they don't exist on one of the platforms we 
needed to support. 
*/ 

char parse_args (int argc, char **argv, char *inputfile, char *outputbase, 
char *darkfile, char *flatfile, float *offset, float *exptime, 
float *numimg, char *cosmicbool, char *readnoise, 
char *cosmicout, float *dccurrent, float *initexp, 
float *shieldcover , char *exptimef ile , float *crrate, 
float *rnvar, float *gain, PIXEL *P) { 
int i ; 

FILE * input; 

/* Intialize various outputs with default values. */ 
*exptime = 1000.; //Exp exptime = 1000s 

*numimg = 64.; //64 readings within the obs time, numbered 0...63. 

strcpy (outputbase , "datasim_"); //Output files named datasim_xxxxxx . f it , 
// where xxxxxx is the integral number of seconds into the observation 
// which the data was taken. 

strcpy (inputf ile , ""); //Default is an empty image (No file), 
strcpy (darkf ile , ""); //Default dark field is an image 

for (i=l; Kargc; i++) { 

if (strcmp(argv[i] , "-h") == 0) { 

printf("%s inputf ile [-t exptime -n numobs -o outputbase -h]\n", 
argv [0] ) ; 
return 1 ; 
} else 

if (strcmp(argv[i] , "-t") == 0) { 

*exptime = string2num(argv [i+1] ) ; 

i+= 1; 
} else 

if (strcmp(argv [i] , "-n") == 0) { 

*numimg = string2num (argv [i+1] ) ; 

i+=l; 
} else 
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if (strcmp(argv[i] , "-o") == 0) { 
strcpy (outputbase , argv[i+l] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-df") == 0) { 
strcpy (darkf ile,argv[i+l] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-ff") == 0) { 
strcpy(f latf ile,argv[i+l] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-de") == 0) { 
*dccurrent = string2num(argv [i+1] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-ie") == 0) { 
*initexp = string2num(argv [i+1] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-sc") == 0) { 

*shieldcover = string2num(argv [i+1] ) ; 
i+=l; 

} else 

if (strcmp(argv[i] , "-co") == 0) { 
*offset = string2num(argv [i+1] ) ; 
} else 

if (strcmp(argv[i] , "-v") == 0) { 
VERBOSE = 1; 
} else 

if (strcmp(argv[i] , "-nc") == 0) { 
*cosmicbool += 1; 
printf ("NO CRs\n") ; 

> 

else 

if (strcmp(argv[i] , "-ncf ") == 0) { 
*cosmicout = 1; 

> 

else 

if (strcmp(argv[i] , "-fs") == 0) { 
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*cosmicbool += 2; 

printf ("Fowler Sampling\n") ; 

} 

else 

if (strcmp(argv[i] , "-expf ") == 0) { 
*cosmicbool += 8; 
strcpy (exptimef ile , argv[++i]); 
printf ("Custom Pattern: °/ s\n" ,argv[i] ) ; 

} 

else 

if (strcmp(argv[i] , "-np") == 0) { 
*cosmicbool += 4; 
printf ("No Poisson Noise\n"); 

} 

else 

if (strcmp(argv[i] , "-craa") == 0) { 
*cosmicbool += 16; 

printf ("Using Cosmic Ray current-leak. \n") ; 

} 

else 

if (strcmp(argv[i] , "-nr") == 0) { 
*readnoise = 1; 
printf ("No read noise\n"); 

} 

else 

if (strcmp(argv[i] , "-err") == 0) { 
*crrate = string2num(argv [++i] ) ; 

} 

else 

if (strcmp(argv[i] , "-r") == 0) { 

*rnvar = string2num(argv [++i] ) ; 
} else 

if (strcmp(argv[i] , "-g") == 0) { 
*gain = string2num(argv [++i] ) ; 
} else 

if (strcmp(argv[i] , "-pix") == 0) { 
char a [100] , *b=argv [++i] ; 

strncpy(a,b,strcspn(b, " , ")) ; 
P->x = string2num(a) ; 
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b += strcspn(b, " , ") +1;; 
strncpy(a,b,strcspn(b, " , ")) ; 
a[strcspn(b,",")] = 0; 
P->y = string2num(a) ; 

b += strcspn(b, " , 
P->z = string2num(b) ; 

} 

else 

strcpy (inputf ile , argv[i] ) ; 

> 

return 0; 

} 



G. string2num.cc 

/ ******* 

string2num. c 

C Code file containing a routine to convert a string to a number (floating 
point) . 

Written by Joel D. Offenberg, Raytheon ITSS 
******/ 

#include <string.h> 
#include <math.h> 

float string2num (char *string) 
{ 

float res=0.0, place=l . ; 
int i , 1 ; 
char dp = 0; 

1 = strlen(string) ; 
for (i=0; i<l; i++) { 



-35- 



switch (string [i] ) { 



CcLSG 


'0' : 


CcLSG 


' 1 ' : 


case 


'2' : 


case 


'3' : 


case 


> 4 > . 


case 


'5' : 


case 


' 6 ' : 


case 


'7' : 


case 


>8> : 


case 


'9' : 


if 


(dp == 1) 


place /= 


10; 



res += (string [i] - '0') * place; 
} 

else { 
res *= 10.0; 

res += (float) string [i] - '0'; 
} 

break; 
case ' . ' : 
dp = 1; 
break; 
case ' e ' : 
case 'E' : 

if (string [i+1] == '-') { 
res /= pow(10 . , string2num(string+i+2) ) ; 

} else { 
float tmp = 0, tmp2; 
tmp = string2num(string+i+2) ; 
tmp2 = pow (10.0, tmp) ; 

res *= pow(10 . , string2num(string+i+2) ) ; 
} 

i=i; 

break; 
default : 
break; 

} 

> 
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return res; 

} 



